clear all
set more off
set mem 10000000
set matsize 10000
version 15

********************************************************************************************
*** RDROBUST reduced-form, Census outcomes: including states with bad/missing shapefiles ***
********************************************************************************************

** Set file paths
do "$path_code/paths.do"

** Set graph scheme
cd "$path/code/analyze"
set scheme fb, perm

****************************************************************** 
****************************************************************** 

{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in RD sample (INCLUDING THE 10 STATES WITH BAD/MISSING SHAPEFILES
gen in_rf_sample = vplan4<11 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_rf_sample==1 
	
	// Grab computer variables from HPCA (didn't get merged into final dataset for some reason
merge m:1 v_id_hpca11 using "$hpca11/hpca11.dta" , keepusing( hpca11_assets_cl_no_int hpca11_assets_cl_int) keep(1 3) nogen
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
drop if st_code==11 // Sikkim, only 5 villages
tab st_code if in_rf_sample==1, gen(STFErf)	
drop STFErf1 STFErf2 STFErf4 // to avoid collinearity (also drop Haryana and Himachal Pradesh FEs, since there are too few observations)

	// Create district FEs (since RD robust doesn't let you pass them through)
gen stdtFE = stdt
tab stdtFE state     if inlist(stdtFE,62,63,64,70,91,175,180,181,202,308,309,320,490,493,494,495,504,505,508)	
replace stdtFE = 999 if inlist(stdtFE,62,63,64,70,91,175,180,181,202,308,309,320,490,493,494,495,504,505,508)
	// one catch-all district FE for districts with so few in-sample villages that they break rdrobust
tab stdtFE if in_rf_sample==1, gen(DTFErf)	
drop DTFErf1 // to avoid collinearity

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

// 	// Create lights-difference variable, to identify crazy outliers
// gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

	// Create population difference outcome variables
gen tot_p11_diff = tot_p11-tot_p
reghdfe tot_p11_diff tot_p if tot_p<1000, a(st_code) residuals(tot_p11_diff2)

	// Create macros for lists of outcome variables 
global vars_demo = "tot_p11 tot_p11_diff tot_p11_diff2 pct_06_11 hpca11_h_size_avg lit_p_11"
global vars_labor_pct = "work_pooled_ag_m_11 work_pooled_ag_f_11 work_pooled_hh_m_11 work_pooled_hh_f_11 work_pooled_ot_m_11 work_pooled_ot_f_11 work_pooled_ag_p_11 work_pooled_ot_p_11 work_main_p_11 work_marg_p_11"
global vars_assets = "hpca11_assets_tel_l_m_b hpca11_assets_tv hpca11_assets_bic hpca11_assets_smm hpca11_assets_none hpca11_assets_rt hpca11_assets_cl_int hpca11_assets_cl_no_int"
global vars_hhold = "hpca11_cook_good hpca11_msl_bad hpca11_mat_f_m hpca11_mat_r_gtbw hpca11_t_hh_d hpca11_msl_elec"
global vars_vd = "vdp_com_d_phone_mob_11 vd_com_d_post_off_11 vd_fin_d_ac_soc_11 vd_wat_d_tubewell_11 pct_irr_11 pct_sown_11"

	// Create index of HPCA variables
gen hpca11_mat_r_not_gtbw = 1-hpca11_mat_r_gtbw	
gen hpca01_mat_r_not_gtbw = 1-hpca01_mat_r_gtbw	
gen hpca11_mat_f_not_m = 1-hpca11_mat_f_m
gen hpca01_mat_f_not_m = 1-hpca01_mat_f_m	
local hpca_index_vars = "hpca11_t_hh_g hpca11_t_hh_r_g hpca11_t_hh_rc_g hpca11_owned hpca11_w_within hpca11_lat_premise hpca11_cook_good " ///
					  + "hpca11_cook_ih_kitchen hpca11_bank hpca11_assets_tel_l_m_b hpca11_assets_tv hpca11_assets_bic hpca11_assets_smm " ///
					  + "hpca11_hh_perm hpca11_mat_r_not_gtbw hpca11_mat_f_not_m " 
foreach v of varlist `hpca_index_vars' {
    local v01 = subinstr("`v'","hpca11","hpca01",1)
	qui sum `v' if tot_p<=1000
	gen Z`v' = (`v' - r(mean))/(r(sd))
	qui sum `v01' if tot_p<=1000
	gen Z`v01' = (`v01' - r(mean))/(r(sd))
}
egen double index_hpca11 = rmean(Zhpca11*)
egen double index_hpca01 = rmean(Zhpca01*)
drop Z*

	// Create index of VD variables
local vd_index_vars = "vd_com_d_phone_ll_11 vd_com_d_post_off_11 vd_edu_n_p_sch_11 vd_edu_n_m_sch_11 vd_edu_n_s_sch_11 " ///
					+ "vd_edu_n_s_s_sch_11 vd_edu_n_coll_11 vd_edu_n_tr_sch_11 vd_fin_d_ac_soc_11 vd_fin_d_bank_11 vd_fin_d_comm_bank_11 " /// 
					+ "vd_fin_d_coop_bank_11 vd_geo_a_irr_11 vd_hea_n_alt_hosp_11 vd_hea_n_fw_cntr_11 vd_hea_n_mcw_cntr_11 vd_hea_n_ph_subcntr_11 " ///
					+ "vd_tra_d_bus_11 vd_wat_d_any_11"
local vd_index_vars11 = "vdp_com_d_phone_mob_11 vdp_tra_d_auto_11 vdp_tra_d_aw_road_11 vdp_tra_d_rick_11 vdp_tra_d_taxi_11 vdp_tra_d_tractor_11 vdp_tra_d_van_11 pct_sown_11"
foreach v of varlist `vd_index_vars' `vd_index_vars11' {
    local v01 = subinstr("`v'","_11","_01",1)
	qui sum `v' if tot_p<=1000
	gen Z`v' = (`v' - r(mean))/(r(sd))
	cap {
		qui sum `v01' if tot_p<=1000
		gen Z`v01' = (`v01' - r(mean))/(r(sd))
	}	
}
egen double index_vd11 = rmean(Zvd*11)
egen double index_vd01 = rmean(Zvd*01)
drop Z*

global vars_index = "work_p_11 index_hpca11 index_vd11"


	// Create variables to store regresson results
gen rf_step = .
gen yvar = ""
gen ifs = ""
gen controls_base = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen nobs_total = .
gen ndist = . 
gen nstate = . 
gen ymean = .
gen ftag = ""
global h_wide = 200


	// Massive loop over RDROBUST outcomes and sensitivities
foreach y in $vars_demo $vars_labor_pct $vars_assets $vars_hhold $vars_vd $vars_index { 
	
	// Prep to store results
	foreach v of varlist rf_step-ftag {
		cap replace `v' = ""
		cap replace `v' = .
	}
	local row = 0

	// Define outcome-specific stuff
	global yvar = "`y'"
	global control = subinstr("$yvar","11","01",1)
	if substr("$yvar",-4,4)=="1101" {
		global control = subinstr("$yvar","1101","01",1)
	}
	do "$path/code/analyze/RDROBUST_rf_outcomes_graphspecs.do"

	//Loop through sensitivities
	foreach rf_step in 1 4 5 11 {

		// Reset RDROBUST defaults
		local ifs = "in_rf_sample==1 & pop_mismatch20==0"
		local controls_base = "" 
		local fe = "STFErf*"
		local kernel = "tri"
		local bwmethod = "mserd"
		local vce = ""
		local poly = 1
		local graphtag = "rf"

		// Define step-specific RDROBUST settings
		if inlist(`rf_step',1) {
			local folder = "RDROBUST plots rf outcomes preferred"
			local ftag = "RF outcomes, preferred, including bad shapfiles"
			local title2 = ""
		}
		if inlist(`rf_step',2) {
			local folder = "RDROBUST plots rf outcomes popmismatch"
			local ftag = "RF outcomes, popmismatch, including bad shapfiles"
			local title2 = "(population mismatches)"
			local ifs = "in_rf_sample==1"
		}
		if inlist(`rf_step',4) {
			local folder = "RDROBUST plots rf outcomes epa kernel"
			local ftag = "RF outcomes, epa kernel, including bad shapfiles"
			local title2 = "(Epanechnikov kernel)"
			local kernel = "epa"
		}
		if inlist(`rf_step',5) {
			local folder = "RDROBUST plots rf outcomes uni kernel"
			local ftag = "RF outcomes, uni kernel, including bad shapfiles"
			local title2 = "(uniform kernel)"
			local kernel = "uni"
		}
		if inlist(`rf_step',7) {
			local folder = "RDROBUST plots rf outcomes no FEs"
			local ftag = "RF outcomes, no FEs, including bad shapfiles"
			local title2 = "(no fixed effects)"
			local fe = ""
		}
		if inlist(`rf_step',8) {
			local folder = "RDROBUST plots rf outcomes district FEs"
			local ftag = "RF outcomes, district FEs, including bad shapfiles"
			local title2 = "(district fixed effects)"
			local fe = "DTFEfs*"
		}
		if inlist(`rf_step',9) {
			local folder = "RDROBUST plots rf outcomes nncluster by district"
			local ftag = "RF outcomes, nncluster by district, including bad shapfiles"
			local title2 = "(nncluster by district)"
			local vce = "vce(nncluster stdt)"
		}
		if inlist(`rf_step',10) {
			local folder = "RDROBUST plots rf outcomes cluster by district"
			local ftag = "RF outcomes, cluster by district, including bad shapfiles"
			local title2 = "(cluster by district)"
			local vce = "vce(cluster stdt)"
		}
		if inlist(`rf_step',11) {
			local folder = "RDROBUST plots rf outcomes CERRD bandwidth"
			local ftag = "RF outcomes, CERRD bandwidth, including bad shapfiles"
			local title2 = "(CERRD bandwidth)"
			local bwmethod = "cerrd"
		}
		
		// Run first-stage regresssion
		di "rdrobust $yvar tot_p if `ifs', c(299.5) covs(`controls_base' $control `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'"
		rdrobust $yvar tot_p if `ifs', c(299.5) covs(`controls_base' $control `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & ${yvar}!=.
		local h_l = e(h_l)
		local h_r = e(h_r)

		// Store results
		local row = `row' + 1
		qui replace rf_step = `rf_step' in `row'
		qui replace yvar = "$yvar" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace controls_base = "`controls_base'" in `row'
		qui replace control = "$control" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui replace nobs_total = e(N_h_l) + e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui unique st_code if temp_in_reg==1
		qui replace nstate = r(unique) in `row'
		qui sum $yvar if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize $yvar for regression sample
		qui reg $yvar `controls_base' $control `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize $yvar for [100,500] bandwidth
		qui reg $yvar `controls_base' $control `fe' if `ifs' & inrange(tot_p,299.5-${h_wide},299.5+${h_wide})
		cap drop y_resid_full
		predict y_resid_full, residuals

		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("$title" "`title2'", color(black) size(large)) ///
				ytitle("$ytitle", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/${yvar}_`graphtag'_inreg_`bins'.pdf", replace
		}

		if !inlist(`rf_step',11) {
			// RD plots for regression sample (constant bandwidth) 
			foreach bins in 20 40 {
				rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-${h_wide},299.5+${h_wide}), ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(${h_wide} ${h_wide}) ///
					graph_options( ///
					title("$title" "`title2'", color(black) size(large)) ///
					ytitle("$ytitle", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(, labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/${yvar}_`graphtag'_wide_`bins'.pdf", replace
			}
		}
	}

	// Save results
	preserve
	keep rf_step-ftag
	dropmiss, obs force
	if "$yvar"!="tot_p11" {
		tempfile reg_results
		save `reg_results'
		clear
		append using "$results/RDROBUST_outcomes_rf_census_badshapefiles.dta" `reg_results'
	}
	duplicates drop
	compress
	save "$results/RDROBUST_outcomes_rf_census_badshapefiles.dta", replace
	restore
	
}



}

****************************************************************** 
****************************************************************** 

